# Load library
library(dplyr)
library(survival)
library(janitor)
library(magrittr)
library(car)
library(ggplot2)
library(tidyverse)
library(broom)
library(MASS)
library(boot)
#print(citation("survival"), bibtex=TRUE)Cox regression modeling of survival after chemotherapy for colon cancer
Data: Survival after chemotherapy for Stage B/C colon cancer [1] [2]
Description
These are data from one of the first successful trials of adjuvant chemotherapy for colon cancer. Levamisole is a low-toxicity compound previously used to treat worm infestations in animals; 5-FU is a moderately toxic (as these things go) chemotherapy agent. There are two records per person, one for recurrence and one for death.
The purpose of this project is to compare survival between the untreated (Obs) group vs those treated with amisole (Lev), or amisole + 5-FU.
Column names:
| id: | id |
| study: | 1 for all patients |
| rx: | Treatment - Obs(ervation), Lev(amisole), Lev(amisole)+5-FU |
| sex: | 1=male |
| age: | in years |
| obstruct: | obstruction of colon by tumour |
| perfor: | perforation of colon |
| adhere: | adherence to nearby organs |
| nodes: | number of lymph nodes with detectable cancer |
| time: | days until event or censoring |
| status: | censoring status |
| differ: | differentiation of tumour (1=well, 2=moderate, 3=poor) |
| extent: | Extent of local spread (1=submucosa, 2=muscle, 3=serosa, 4=contiguous structures) |
| surg: | time from surgery to registration (0=short, 1=long) |
| node4: | more than 4 positive lymph nodes |
| etype: | event type: 1=recurrence,2=death |
#Load data
colon <- as_tibble(colon)
head(colon)# A tibble: 6 × 16
id study rx sex age obstruct perfor adhere nodes status differ
<dbl> <dbl> <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 1 1 Lev+5FU 1 43 0 0 0 5 1 2
2 1 1 Lev+5FU 1 43 0 0 0 5 1 2
3 2 1 Lev+5FU 1 63 0 0 0 1 0 2
4 2 1 Lev+5FU 1 63 0 0 0 1 0 2
5 3 1 Obs 0 71 0 0 1 7 1 2
6 3 1 Obs 0 71 0 0 1 7 1 2
# ℹ 5 more variables: extent <dbl>, surg <dbl>, node4 <dbl>, time <dbl>,
# etype <dbl>Since the current analysis is focused on survival, filter data to death as the event type. This will create a data table with one row per individual.
colon_surv <- colon%>%filter(etype == 2) Identify participants who had recurrence. Identify those not censored for recurrence event. Filter event type = 1 (recurrence), status = 1.
recurrence <- colon%>%filter(etype == 1 & status == 1)%>%dplyr::select(id)
recurrence <- recurrence%>%mutate(recurrence = 1) # list of patients with recurrence
colon_surv <- colon_surv%>%merge(recurrence, by = "id", all.x = TRUE)
colon_surv$recurrence[is.na(colon_surv$recurrence)] <- 0I. Exploratory data analysis
Check missing values
na_counts <- sapply(colon_surv, function(x)sum(is.na(x)))
na_counts id study rx sex age obstruct perfor
0 0 0 0 0 0 0
adhere nodes status differ extent surg node4
0 18 0 23 0 0 0
time etype recurrence
0 0 0 # replace NAs with mode
table(colon_surv$differ)
1 2 3
93 663 150 mode(colon_surv$differ)[1] "numeric"median(colon_surv$nodes, na.rm= TRUE)[1] 2colon_surv$differ <- if_else(is.na(colon_surv$differ), 2,colon_surv$differ)
colon_surv$nodes <- if_else(is.na(colon_surv$nodes), 2,colon_surv$nodes)Insight: only nodes and differ columns have NA values. Replacing the 23 NAs in differ column with mode, and replace NAs in nodes with median.
Evaluate continuous variables
# age
hist(colon_surv$age)
hist(colon_surv$nodes)
hist(colon_surv$time)
Insight: Age is normally distribute. Number of nodes is skewed to the right. Time is fairly normally distributed with most the individuals had event time between 500-3000 days.
Evaluate nodes column to investigate outliers
t <- colon_surv%>%filter(node4 ==1) # samples with more than positive lymph nodes
hist(t$nodes) 
Insight: samples with greater than 4 lymph nodes have less than 5 count in nodes column, so the two columns are not consistent. Therefore, nodes column will not be used for further analysis.
Evaluate categorical variables
summary_table <- colon_surv%>%summarise(count =n(),
male = sum(sex),
median_age = median(age),
ct_perforation = sum(perfor),
ct_adherence_nerby_organ = sum(adhere))
summary_table count male median_age ct_perforation ct_adherence_nerby_organ
1 929 484 61 27 135Insight: Total number of participants: 929. About half of the participants are male and about half were censored, while the other half died.
colon_surv <- colon_surv%>%mutate(differentiation = case_when(differ == 1 ~ "well",
differ == 2 ~ "moderate",
differ == 3 ~ "poor"),
local_spread = case_when(extent == 1 ~ "submucosa",
extent == 2 ~ "muscle",
extent == 3 ~ "serosa",
extent == 4 ~ "contiguous"),
surg_to_reg_time = case_when(surg == 0~ "short",
surg == 1 ~ "long"))Frequency tables for categorical variables
# frequency tables for categorical variables
# Tumor differentiation
colon_surv %>%
tabyl(differentiation, rx) %>%
adorn_percentages("col") %>%
adorn_pct_formatting(digits = 1) %>%
adorn_ns() differentiation Obs Lev Lev+5FU
moderate 74.9% (236) 73.9% (229) 72.7% (221)
poor 16.5% (52) 14.2% (44) 17.8% (54)
well 8.6% (27) 11.9% (37) 9.5% (29)# extent of local spread
colon_surv %>%
tabyl(local_spread, rx) %>%
adorn_percentages("col") %>%
adorn_pct_formatting(digits = 1) %>%
adorn_ns() local_spread Obs Lev Lev+5FU
contiguous 6.3% (20) 3.9% (12) 3.6% (11)
muscle 12.1% (38) 11.6% (36) 10.5% (32)
serosa 79.0% (249) 83.5% (259) 82.6% (251)
submucosa 2.5% (8) 1.0% (3) 3.3% (10)# colum obstruction
colon_surv %>%
tabyl(obstruct, rx) %>%
adorn_percentages("col") %>%
adorn_pct_formatting(digits = 1) %>%
adorn_ns() obstruct Obs Lev Lev+5FU
0 80.0% (252) 79.7% (247) 82.2% (250)
1 20.0% (63) 20.3% (63) 17.8% (54)# colon perforation
colon_surv %>%
tabyl(perfor, rx) %>%
adorn_percentages("col") %>%
adorn_pct_formatting(digits = 1) %>%
adorn_ns() perfor Obs Lev Lev+5FU
0 97.1% (306) 96.8% (300) 97.4% (296)
1 2.9% (9) 3.2% (10) 2.6% (8)# Adherance to nearby organs
colon_surv %>%
tabyl(adhere, rx) %>%
adorn_percentages("col") %>%
adorn_pct_formatting(digits = 1) %>%
adorn_ns() adhere Obs Lev Lev+5FU
0 85.1% (268) 84.2% (261) 87.2% (265)
1 14.9% (47) 15.8% (49) 12.8% (39)# extent of local tumor spread
colon_surv %>%
tabyl(local_spread, rx) %>%
adorn_percentages("col") %>%
adorn_pct_formatting(digits = 1) %>%
adorn_ns() local_spread Obs Lev Lev+5FU
contiguous 6.3% (20) 3.9% (12) 3.6% (11)
muscle 12.1% (38) 11.6% (36) 10.5% (32)
serosa 79.0% (249) 83.5% (259) 82.6% (251)
submucosa 2.5% (8) 1.0% (3) 3.3% (10)# More than 4 lymph nodes with cancer
colon_surv %>%
tabyl(node4, rx) %>%
adorn_percentages("col") %>%
adorn_pct_formatting(digits = 1) %>%
adorn_ns() node4 Obs Lev Lev+5FU
0 72.4% (228) 71.3% (221) 74.0% (225)
1 27.6% (87) 28.7% (89) 26.0% (79)# Recurrence
colon_surv %>%
tabyl(recurrence, rx) %>%
adorn_percentages("col") %>%
adorn_pct_formatting(digits = 1) %>%
adorn_ns() recurrence Obs Lev Lev+5FU
0 43.8% (138) 44.5% (138) 60.9% (185)
1 56.2% (177) 55.5% (172) 39.1% (119)colon_surv %>%
tabyl(recurrence, status) %>%
adorn_percentages("col") %>%
adorn_pct_formatting(digits = 1) %>%
adorn_ns() recurrence 0 1
0 88.7% (423) 8.4% (38)
1 11.3% (54) 91.6% (414)# time from surgery to registration
colon_surv %>%
tabyl(surg, rx) %>%
adorn_percentages("col") %>%
adorn_pct_formatting(digits = 1) %>%
adorn_ns() surg Obs Lev Lev+5FU
0 71.1% (224) 74.2% (230) 75.0% (228)
1 28.9% (91) 25.8% (80) 25.0% (76)
Summary statistics grouped by treatment
summary_table <- colon_surv%>%group_by(rx)%>%summarise(count =n(),
male = sum(sex),
median_age = median(age),
ct_perforation = sum(perfor),
ct_adherence_nerby_organ = sum(adhere),
perc_male = (male/count)*100,
iqr_age = IQR(age))
summary_table# A tibble: 3 × 8
rx count male median_age ct_perforation ct_adherence_nerby_o…¹ perc_male
<fct> <int> <dbl> <dbl> <dbl> <dbl> <dbl>
1 Obs 315 166 60 9 47 52.7
2 Lev 310 177 61 10 49 57.1
3 Lev+5FU 304 141 62 8 39 46.4
# ℹ abbreviated name: ¹ct_adherence_nerby_organ
# ℹ 1 more variable: iqr_age <dbl>g <- colon_surv%>%filter(rx == "Lev+5FU")
summary(g$age) Min. 1st Qu. Median Mean 3rd Qu. Max.
26.0 52.0 62.0 59.7 70.0 81.0 Insight: Each treatment group had about 300 participants. Median age, number of participants with perforation and adherence are similar between the three groups.
II. Table 1: Description of the study population
| Observation (%) | Amisole (%) | Amisole + 5-FU (%) | ||
|---|---|---|---|---|
| N=315 | N=310 | N=304 | ||
| Demographics | ||||
| Male | 166 (52.3) | 177 (57.1) | 141 | |
| Median age (years) [IQR] | 60 [53,68] | 61 [53,69] | 61 [52,70] | |
| Cancer characteristics | ||||
| Colon obstruction | 63 (20.0) | 63 (20.3) | 54 (17.8) | |
| Colon perforation | 9 (2.9) | 10 (3.2) | 8 (2.6) | |
| Adherence to nearby organs | 47 (14.9) | 49 (15.8) | 39 (12.8) | |
| Differentiation of tumor | ||||
| Well | 27 (8.6) | 37 (11.9) | 29 (9.5) | |
| Moderate | 236 (74.9) | 229 (73.9) | 221 (72.7) | |
| Poor | 52 (16.5) | 44 (14.2) | 54 (17.8) | |
| Extent of local spread | ||||
| Contiguous | 20 (6.3) | 12 (3.9) | 11 (3.6) | |
| Muscle | 38 (12.1) | 36 (11.6) | 32 (10.5) | |
| Serosa | 249 (79.0) | 259 (83.5) | 251 (82.6) | |
| Submucosa | 8 (2.5) | 3 (1.0) | 10 (3.3) | |
| More than 4 lymph nodes with cancer | Yes | 87 (27.6) | 89 (28.7) | 79 (26.0) |
| Recurrence (%) | Yes | 177 (56.2) | 172 (55.5) | 119 (39.1) |
| Short time from surgery to registration (%) | Yes | 91 (28.9) | 80 (25.8) | 76 (25.0) |
III. Methods
The Cox proportional hazards model was used to model the relationship between survival time and different lung cancer treatments. In particlular the survival time will be compared between the untreated group (observation) vs. those treated with amisole (Lev), or amisole + 5-FU. The Cox regression model was chosen for this study because it is useful for studying association between survival time of patients and predictors and allows estimating the relative risk or hazard ratios due to the covariates, i.e., treatment status. The time (in days) until event, i.e, death, will be modeled as a function of treatment and other variables, including age, sex, and various tumor characteristics. Significant predictors were included in the final model.
Statistical analysis
The R statistical software version 4.3.2 [3] was used for all analysis. The Survival package was used to construct the Cox regression model [2] [1].
Cox regression model is based on the hazard function \(h_x(t)\) with covariates at time t given by:
\(h_x(t)=h_0(t)\exp(\beta_1x_1 +\beta_2x_2 + \dots + \beta_p x_p)\)
Where:
\(h_x(t)\) is the hazard function
\(h_0(t)\) is the baseline hazard function
\(\beta_1x_1 + \beta_2x_2 + \dots +\beta_p x_p\) represent the linear combination of covariates and their coefficient
The hazards for the observation vs. amisole (Lev), or amisole + 5-FU group with covariate values x1 and x2 are given by: \(hx_1(t)=h_0(t)\exp(\beta_1x_1)\) and \(hx_2(t)=h_0(t)\exp(\beta_2x_2)\), respectively
The hazard ratio is expressed as: HR = \(hx_2(t)\) / \(hx_1(t)\) = \(\exp[\beta(x_2-x_1)]\)
The Schoenfeld residual plot was constructed to test Cox proportional hazards assumption. When the proportional hazards assumpiton was not met for any of the covariates, stratification approach was explored. The Survminer [4] package was used to plot the Kaplan-Meier curve to visualize the survival probability over time for each treatment group.
Multicolinearity was tested using Variant Inflation Factor (VIF) calculated using MASS package [5].
The R MASS package was used for Stepwise model selection, using “both” forward and backward variable selection [5]. For Stepwise selection, stepAIC() function uses AIC (Akaike Information Criterion) as the measure to add or remove predictors from the model. Model performance was evaluated using 100-fold cross-validation using the boot package [6] [7].
IV. Analysis: Cox regression model
Plot survival times
# Create new incremental count id
colon_surv$idcount <- c(1:length(colon_surv$id))
# Order by survival time and create an order variable:
colon_surv <- colon_surv[order(-colon_surv$time, colon_surv$status),]
colon_surv$order <- c(1:length(colon_surv$idcount))
ggplot(data=colon_surv, aes(x=time, y=order)) +
geom_rect(xmin=23,xmax=colon_surv$time,ymin=colon_surv$order,ymax=colon_surv$order+1, colour="lightgray") +
geom_rect(xmin=colon_surv$time-2,xmax=colon_surv$time,ymin=colon_surv$order,ymax=colon_surv$order+1,
fill=factor(colon_surv$status+1)) +
geom_vline(xintercept= 1976,linetype="solid") +
scale_x_continuous(breaks=seq(20,3330,650)) +
geom_text(aes(2600, 750, label="Median Survival Time")) +
xlab("Survival Time (Days)") + ylab("Participants (ordered by survival time)") +
ggtitle("Survival Times for Participant") +
theme_classic() +
theme(legend.position="none",
panel.grid.major=element_blank(),
panel.grid.minor=element_blank(),
panel.background=element_blank(),
axis.line.x = element_line(color = "black"),
axis.line.y = element_line(color = "black"))
Plot survival curve stratified by treatment group
library(survminer)
library(survival)
# Estimate the median survival time among the three groups
survfit(Surv(time,status) ~ rx, data = colon_surv)Call: survfit(formula = Surv(time, status) ~ rx, data = colon_surv)
n events median 0.95LCL 0.95UCL
rx=Obs 315 168 2083 1656 2789
rx=Lev 310 161 2152 1540 NA
rx=Lev+5FU 304 123 NA 2725 NA# count the number of events after 2080 days, which is the median survival time among the observation group
tt <- colon_surv%>%filter(time > 2083)%>% group_by(rx)%>%summarise(ct = n(),
death = sum(status))
# Plot survival curve
fit <- survfit(Surv(time,status) ~ rx, data = colon_surv)
ggsurvplot(fit, data=colon_surv, risk.table = TRUE)
# Estimate the probability of surviving beyond 3000 days
summary(survfit(Surv(time, status) ~ rx, data = colon_surv), times = 3000)Call: survfit(formula = Surv(time, status) ~ rx, data = colon_surv)
rx=Obs
time n.risk n.event survival std.err lower 95% CI
3.00e+03 6.00e+00 1.68e+02 4.08e-01 3.97e-02 3.37e-01
upper 95% CI
4.94e-01
rx=Lev
time n.risk n.event survival std.err lower 95% CI
3.00e+03 4.00e+00 1.61e+02 3.92e-01 5.63e-02 2.96e-01
upper 95% CI
5.20e-01
rx=Lev+5FU
time n.risk n.event survival std.err lower 95% CI
3.00e+03 7.00e+00 1.23e+02 5.61e-01 3.43e-02 4.97e-01
upper 95% CI
6.32e-01 # compare significant diffeerence is survival times between the three groups
survdiff(Surv(time, status)~ rx, data = colon_surv)Call:
survdiff(formula = Surv(time, status) ~ rx, data = colon_surv)
N Observed Expected (O-E)^2/E (O-E)^2/V
rx=Obs 315 168 148 2.58 3.85
rx=Lev 310 161 146 1.52 2.25
rx=Lev+5FU 304 123 157 7.55 11.62
Chisq= 11.7 on 2 degrees of freedom, p= 0.003 Insight: Based on the survival curve, the mediant survival time for the observation group is 2083 days. However, the median survival of Lev and Lev+5Fu group cannot be estimated, because there are too few events after 2083 days, which is the median survival time in the observation group.
The time for 50% survival probability of the group treated with Lev+5Fu is over 3000 days while the survival time for the observation and Lev group is around 2080 days. The probability of surviving to 3000 days among the Lev+5FU group is 56% (95% CI: 50-63), compared to 41% among the observation group.
The survival time is significantly different (P=0.003) between the three groups.
Fit base Model
m0 <- coxph(Surv(time, status) ~ 1, data = colon_surv)
summary(m0)Call: coxph(formula = Surv(time, status) ~ 1, data = colon_surv)
Null model
log likelihood= -2930.192
n= 929 Univariate analysis
Test significance of treatment as a predictor
#| echo: true
#| message: false
#| warning: false
# Univariate analysis
m1 <- coxph(Surv(time, status) ~ rx, data = colon_surv)
summary(m1)Call:
coxph(formula = Surv(time, status) ~ rx, data = colon_surv)
n= 929, number of events= 452
coef exp(coef) se(coef) z Pr(>|z|)
rxLev -0.02664 0.97371 0.11030 -0.241 0.80917
rxLev+5FU -0.37171 0.68955 0.11875 -3.130 0.00175 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
exp(coef) exp(-coef) lower .95 upper .95
rxLev 0.9737 1.027 0.7844 1.2087
rxLev+5FU 0.6896 1.450 0.5464 0.8703
Concordance= 0.536 (se = 0.013 )
Likelihood ratio test= 12.15 on 2 df, p=0.002
Wald test = 11.56 on 2 df, p=0.003
Score (logrank) test = 11.68 on 2 df, p=0.003Insight: the coefficient of Lev is not significant, suggesting that there is no evidence that this treatment affects survival time compared to observation. however Lev+5Fu is significant (p=0.00175), indicating that the treatment Lev +5Fu has a statistically significant effect on survival time compared to the reference group. The negative sign indicates that this treatment group has a lower hazard and likely a longer survival time.
The hazard ratio for Lex+5FU (0.690), indicating the risk of death is about 31% lower compared to the observation group.
The p-values indicate that the model is significant.
Test the Cox proportional hazard assumption
cox.zph(m1) chisq df p
rx 1.48 2 0.48
GLOBAL 1.48 2 0.48zph_test <- cox.zph(m1)
print(zph_test) chisq df p
rx 1.48 2 0.48
GLOBAL 1.48 2 0.48# plot the Schoenfeld residuals
plot(zph_test)
Insight: The Schoenfeld residal plot shows that the residuals are scattered randomly and the smooth trend line is horizontal near 0. This suggests that the hazard ratio for rx (treatment status) is constant over time and the proportional hazard assumption is met. The global p-value is >0.05, indicating that the the assumption is met.
Multivariate analysis
Include all variables to determine which predictors are significant.
# Subset data for modeling
df <- colon_surv%>%dplyr::select(!c(id,study,etype,differ,recurrence, extent,surg_to_reg_time, idcount, order, nodes))
# multivariant analysis
m2 <- coxph(Surv(time, status) ~ rx+ age + sex + perfor + adhere + surg + obstruct + differentiation + node4+
local_spread, data = df)
summary(m2)Call:
coxph(formula = Surv(time, status) ~ rx + age + sex + perfor +
adhere + surg + obstruct + differentiation + node4 + local_spread,
data = df)
n= 929, number of events= 452
coef exp(coef) se(coef) z Pr(>|z|)
rxLev -0.0160445 0.9840835 0.1115300 -0.144 0.885612
rxLev+5FU -0.3742695 0.6877915 0.1196533 -3.128 0.001760 **
age 0.0070459 1.0070708 0.0040596 1.736 0.082633 .
sex 0.0412783 1.0421421 0.0952633 0.433 0.664791
perfor 0.0005371 1.0005373 0.2695697 0.002 0.998410
adhere 0.1689086 1.1840119 0.1295242 1.304 0.192210
surg 0.2362177 1.2664500 0.1032961 2.287 0.022207 *
obstruct 0.2865577 1.3318350 0.1173103 2.443 0.014576 *
differentiationpoor 0.3579964 1.4304604 0.1222148 2.929 0.003398 **
differentiationwell 0.0812878 1.0846830 0.1662752 0.489 0.624930
node4 0.9353140 2.5480134 0.0988712 9.460 < 2e-16 ***
local_spreadmuscle -0.9503933 0.3865890 0.2554577 -3.720 0.000199 ***
local_spreadserosa -0.4526719 0.6359268 0.1986566 -2.279 0.022687 *
local_spreadsubmucosa -1.2432513 0.2884449 0.5396100 -2.304 0.021224 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
exp(coef) exp(-coef) lower .95 upper .95
rxLev 0.9841 1.0162 0.7909 1.2245
rxLev+5FU 0.6878 1.4539 0.5440 0.8696
age 1.0071 0.9930 0.9991 1.0151
sex 1.0421 0.9596 0.8646 1.2561
perfor 1.0005 0.9995 0.5899 1.6970
adhere 1.1840 0.8446 0.9186 1.5262
surg 1.2664 0.7896 1.0343 1.5507
obstruct 1.3318 0.7508 1.0583 1.6761
differentiationpoor 1.4305 0.6991 1.1258 1.8176
differentiationwell 1.0847 0.9219 0.7830 1.5026
node4 2.5480 0.3925 2.0991 3.0929
local_spreadmuscle 0.3866 2.5867 0.2343 0.6378
local_spreadserosa 0.6359 1.5725 0.4308 0.9387
local_spreadsubmucosa 0.2884 3.4669 0.1002 0.8306
Concordance= 0.674 (se = 0.012 )
Likelihood ratio test= 147 on 14 df, p=<2e-16
Wald test = 150.3 on 14 df, p=<2e-16
Score (logrank) test = 161.3 on 14 df, p=<2e-16# Determine significant predictors
anova(m2)Analysis of Deviance Table
Cox model: response is Surv(time, status)
Terms added sequentially (first to last)
loglik Chisq Df Pr(>|Chi|)
NULL -2930.2
rx -2924.1 12.1478 2 0.0023022 **
age -2923.9 0.3443 1 0.5573434
sex -2923.9 0.0000 1 0.9964141
perfor -2923.8 0.3345 1 0.5630435
adhere -2921.3 5.0032 1 0.0253001 *
surg -2919.1 4.3247 1 0.0375630 *
obstruct -2916.9 4.4808 1 0.0342771 *
differentiation -2909.8 14.0968 2 0.0008688 ***
node4 -2865.5 88.5800 1 < 2.2e-16 ***
local_spread -2856.7 17.6698 3 0.0005145 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Insight: When all variables are included in the model, the anova test indicates that rx, adhere, surg, obstruct, differentiation, node4 and local spread are significant predictors.
The concordance of the multivariable model, 0.674, is higher than the univariate model (m1, concordance =0.53), suggesting that the multivariate model is a better fit model.
Calculate Variance Inflation Factor (VIF) to assess multicollinearity among predictors
vif <- vif(m2)
print(vif) GVIF Df GVIF^(1/(2*Df))
rx 1.035127 2 1.008668
age 1.042883 1 1.021217
sex 1.022620 1 1.011247
perfor 1.053420 1 1.026363
adhere 1.092721 1 1.045333
surg 1.013876 1 1.006914
obstruct 1.055334 1 1.027295
differentiation 1.062714 2 1.015323
node4 1.046133 1 1.022806
local_spread 1.091169 3 1.014648Insight: None of the variables have VIF values above 5, therefore there is no multicollinearity
Evaluate significance of predictors. Model survival while including different cancer characteristics as predictors separately to identify significance predictors.
# model including all variables
m2 <- coxph(Surv(time, status) ~ rx+ age + sex + perfor + adhere + surg + obstruct + differentiation + node4+
local_spread, data = df)
# Treatment
m2a <- coxph(Surv(time, status) ~ rx, data = colon_surv) # significant
summary(m2a)Call:
coxph(formula = Surv(time, status) ~ rx, data = colon_surv)
n= 929, number of events= 452
coef exp(coef) se(coef) z Pr(>|z|)
rxLev -0.02664 0.97371 0.11030 -0.241 0.80917
rxLev+5FU -0.37171 0.68955 0.11875 -3.130 0.00175 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
exp(coef) exp(-coef) lower .95 upper .95
rxLev 0.9737 1.027 0.7844 1.2087
rxLev+5FU 0.6896 1.450 0.5464 0.8703
Concordance= 0.536 (se = 0.013 )
Likelihood ratio test= 12.15 on 2 df, p=0.002
Wald test = 11.56 on 2 df, p=0.003
Score (logrank) test = 11.68 on 2 df, p=0.003# Demographics
m2b <- coxph(Surv(time, status) ~ age + sex, data = colon_surv) # not significant
summary(m2b)Call:
coxph(formula = Surv(time, status) ~ age + sex, data = colon_surv)
n= 929, number of events= 452
coef exp(coef) se(coef) z Pr(>|z|)
age 0.001951 1.001952 0.004031 0.484 0.628
sex 0.012955 1.013040 0.094205 0.138 0.891
exp(coef) exp(-coef) lower .95 upper .95
age 1.002 0.9981 0.9941 1.010
sex 1.013 0.9871 0.8422 1.218
Concordance= 0.511 (se = 0.014 )
Likelihood ratio test= 0.26 on 2 df, p=0.9
Wald test = 0.25 on 2 df, p=0.9
Score (logrank) test = 0.25 on 2 df, p=0.9# cancer characteristics
m2c <- coxph(Surv(time, status) ~ perfor + adhere + obstruct, data = colon_surv) # adhere and obstruct are significant
summary(m2c)Call:
coxph(formula = Surv(time, status) ~ perfor + adhere + obstruct,
data = colon_surv)
n= 929, number of events= 452
coef exp(coef) se(coef) z Pr(>|z|)
perfor -0.03415 0.96643 0.26932 -0.127 0.8991
adhere 0.31011 1.36358 0.12572 2.467 0.0136 *
obstruct 0.25794 1.29426 0.11547 2.234 0.0255 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
exp(coef) exp(-coef) lower .95 upper .95
perfor 0.9664 1.0347 0.5701 1.638
adhere 1.3636 0.7334 1.0658 1.745
obstruct 1.2943 0.7726 1.0321 1.623
Concordance= 0.536 (se = 0.012 )
Likelihood ratio test= 10.82 on 3 df, p=0.01
Wald test = 11.51 on 3 df, p=0.009
Score (logrank) test = 11.6 on 3 df, p=0.009# Differentiation of tumor
m2d <- coxph(Surv(time, status) ~ differentiation, data = colon_surv) # significant
summary(m2d)Call:
coxph(formula = Surv(time, status) ~ differentiation, data = colon_surv)
n= 929, number of events= 452
coef exp(coef) se(coef) z Pr(>|z|)
differentiationpoor 0.48265 1.62036 0.12040 4.009 6.1e-05 ***
differentiationwell -0.05027 0.95097 0.16408 -0.306 0.759
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
exp(coef) exp(-coef) lower .95 upper .95
differentiationpoor 1.620 0.6171 1.2798 2.052
differentiationwell 0.951 1.0516 0.6895 1.312
Concordance= 0.544 (se = 0.012 )
Likelihood ratio test= 15.34 on 2 df, p=5e-04
Wald test = 16.97 on 2 df, p=2e-04
Score (logrank) test = 17.31 on 2 df, p=2e-04# Extent of local spread
m2f <- coxph(Surv(time, status) ~ local_spread, data = colon_surv) # significant
summary(m2f)Call:
coxph(formula = Surv(time, status) ~ local_spread, data = colon_surv)
n= 929, number of events= 452
coef exp(coef) se(coef) z Pr(>|z|)
local_spreadmuscle -1.0892 0.3365 0.2498 -4.361 1.29e-05 ***
local_spreadserosa -0.5043 0.6039 0.1927 -2.617 0.00888 **
local_spreadsubmucosa -1.7283 0.1776 0.5336 -3.239 0.00120 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
exp(coef) exp(-coef) lower .95 upper .95
local_spreadmuscle 0.3365 2.972 0.2062 0.5490
local_spreadserosa 0.6039 1.656 0.4139 0.8811
local_spreadsubmucosa 0.1776 5.631 0.0624 0.5054
Concordance= 0.552 (se = 0.009 )
Likelihood ratio test= 29.21 on 3 df, p=2e-06
Wald test = 25.35 on 3 df, p=1e-05
Score (logrank) test = 26.94 on 3 df, p=6e-06# Recurrence
# m2g <- coxph(Surv(time, status) ~ recurrence, data = colon_surv) # significant
# summary(m2g)
# short time from surgery to registration
m2h <- coxph(Surv(time, status) ~ surg, data = colon_surv) # significant
summary(m2h)Call:
coxph(formula = Surv(time, status) ~ surg, data = colon_surv)
n= 929, number of events= 452
coef exp(coef) se(coef) z Pr(>|z|)
surg 0.2333 1.2627 0.1026 2.274 0.0229 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
exp(coef) exp(-coef) lower .95 upper .95
surg 1.263 0.7919 1.033 1.544
Concordance= 0.523 (se = 0.011 )
Likelihood ratio test= 5.01 on 1 df, p=0.03
Wald test = 5.17 on 1 df, p=0.02
Score (logrank) test = 5.2 on 1 df, p=0.02# include significant predictors in final model
m2i <- coxph(Surv(time, status) ~ rx + adhere + surg + obstruct + differentiation
+ local_spread, data = colon_surv)
summary(m2i)Call:
coxph(formula = Surv(time, status) ~ rx + adhere + surg + obstruct +
differentiation + local_spread, data = colon_surv)
n= 929, number of events= 452
coef exp(coef) se(coef) z Pr(>|z|)
rxLev -0.0047151 0.9952960 0.1110663 -0.042 0.966137
rxLev+5FU -0.3588384 0.6984872 0.1191372 -3.012 0.002596 **
adhere 0.1499547 1.1617816 0.1285420 1.167 0.243380
surg 0.2103485 1.2341081 0.1031988 2.038 0.041521 *
obstruct 0.2000040 1.2214077 0.1148937 1.741 0.081723 .
differentiationpoor 0.4507752 1.5695283 0.1217364 3.703 0.000213 ***
differentiationwell -0.0003827 0.9996174 0.1651661 -0.002 0.998151
local_spreadmuscle -0.9927505 0.3705561 0.2551148 -3.891 9.97e-05 ***
local_spreadserosa -0.4291507 0.6510618 0.1987757 -2.159 0.030853 *
local_spreadsubmucosa -1.5035698 0.2223351 0.5385347 -2.792 0.005239 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
exp(coef) exp(-coef) lower .95 upper .95
rxLev 0.9953 1.0047 0.80059 1.2373
rxLev+5FU 0.6985 1.4317 0.55303 0.8822
adhere 1.1618 0.8607 0.90304 1.4947
surg 1.2341 0.8103 1.00812 1.5108
obstruct 1.2214 0.8187 0.97513 1.5299
differentiationpoor 1.5695 0.6371 1.23637 1.9925
differentiationwell 0.9996 1.0004 0.72318 1.3817
local_spreadmuscle 0.3706 2.6986 0.22475 0.6110
local_spreadserosa 0.6511 1.5360 0.44099 0.9612
local_spreadsubmucosa 0.2223 4.4977 0.07738 0.6389
Concordance= 0.618 (se = 0.013 )
Likelihood ratio test= 63.62 on 10 df, p=7e-10
Wald test = 60.86 on 10 df, p=2e-09
Score (logrank) test = 63.16 on 10 df, p=9e-10Insight: rx, adhere, surg, obstruct, differentiation and local_spread are significant predictors. However, the model concordance is low (~0.5) when ecah were included separately. Model m2i has the highest conncorance.
Perform Stepwise variable selection:
Use the MASS package stepAIC() function for stepwise selection by using AIC (Akaike Information Criterion) as the measure to add or remove predictors from the model.
library(MASS) # for stepwise regression
# model including all variables
m2 <- coxph(Surv(time, status) ~ rx+ age + sex + perfor + adhere + surg + obstruct + differentiation + node4+
local_spread, data = df)
# stepwise selection
stepwise_model <- stepAIC(m2, direction = "both")Start: AIC=5741.4
Surv(time, status) ~ rx + age + sex + perfor + adhere + surg +
obstruct + differentiation + node4 + local_spread
Df AIC
- perfor 1 5739.4
- sex 1 5739.6
- adhere 1 5741.0
<none> 5741.4
- age 1 5742.5
- surg 1 5744.5
- obstruct 1 5745.1
- differentiation 2 5745.4
- rx 2 5749.9
- local_spread 3 5753.1
- node4 1 5822.0
Step: AIC=5739.4
Surv(time, status) ~ rx + age + sex + adhere + surg + obstruct +
differentiation + node4 + local_spread
Df AIC
- sex 1 5737.6
- adhere 1 5739.1
<none> 5739.4
- age 1 5740.5
+ perfor 1 5741.4
- surg 1 5742.5
- obstruct 1 5743.2
- differentiation 2 5743.5
- rx 2 5747.9
- local_spread 3 5751.1
- node4 1 5820.1
Step: AIC=5737.59
Surv(time, status) ~ rx + age + adhere + surg + obstruct + differentiation +
node4 + local_spread
Df AIC
- adhere 1 5737.3
<none> 5737.6
- age 1 5738.6
+ sex 1 5739.4
+ perfor 1 5739.6
- surg 1 5740.7
- obstruct 1 5741.2
- differentiation 2 5741.7
- rx 2 5746.2
- local_spread 3 5749.3
- node4 1 5818.2
Step: AIC=5737.26
Surv(time, status) ~ rx + age + surg + obstruct + differentiation +
node4 + local_spread
Df AIC
<none> 5737.3
+ adhere 1 5737.6
- age 1 5738.6
+ sex 1 5739.1
+ perfor 1 5739.2
- surg 1 5740.7
- obstruct 1 5740.9
- differentiation 2 5742.1
- rx 2 5746.0
- local_spread 3 5750.4
- node4 1 5817.4summary(stepwise_model)Call:
coxph(formula = Surv(time, status) ~ rx + age + surg + obstruct +
differentiation + node4 + local_spread, data = df)
n= 929, number of events= 452
coef exp(coef) se(coef) z Pr(>|z|)
rxLev -0.010789 0.989269 0.111379 -0.097 0.92283
rxLev+5FU -0.375746 0.686776 0.119575 -3.142 0.00168 **
age 0.007366 1.007393 0.004051 1.818 0.06901 .
surg 0.243871 1.276179 0.103202 2.363 0.01813 *
obstruct 0.283165 1.327324 0.116138 2.438 0.01476 *
differentiationpoor 0.373783 1.453222 0.121599 3.074 0.00211 **
differentiationwell 0.069022 1.071459 0.165690 0.417 0.67699
node4 0.929854 2.534140 0.098507 9.440 < 2e-16 ***
local_spreadmuscle -0.995556 0.369518 0.252106 -3.949 7.85e-05 ***
local_spreadserosa -0.501169 0.605822 0.194154 -2.581 0.00984 **
local_spreadsubmucosa -1.322018 0.266597 0.536289 -2.465 0.01370 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
exp(coef) exp(-coef) lower .95 upper .95
rxLev 0.9893 1.0108 0.79526 1.2306
rxLev+5FU 0.6868 1.4561 0.54329 0.8682
age 1.0074 0.9927 0.99943 1.0154
surg 1.2762 0.7836 1.04248 1.5623
obstruct 1.3273 0.7534 1.05711 1.6666
differentiationpoor 1.4532 0.6881 1.14506 1.8443
differentiationwell 1.0715 0.9333 0.77436 1.4826
node4 2.5341 0.3946 2.08921 3.0738
local_spreadmuscle 0.3695 2.7062 0.22545 0.6057
local_spreadserosa 0.6058 1.6506 0.41408 0.8864
local_spreadsubmucosa 0.2666 3.7510 0.09319 0.7627
Concordance= 0.672 (se = 0.012 )
Likelihood ratio test= 145.1 on 11 df, p=<2e-16
Wald test = 149.1 on 11 df, p=<2e-16
Score (logrank) test = 159.7 on 11 df, p=<2e-16# Multivariate model including variables selected based on stepwise variable selection. The same variables were significant based on anova test of the model that included all variables.
m3 <- coxph(Surv(time, status) ~ rx + age + surg + obstruct +
differentiation + node4 + local_spread, data = df)
summary(m3)Call:
coxph(formula = Surv(time, status) ~ rx + age + surg + obstruct +
differentiation + node4 + local_spread, data = df)
n= 929, number of events= 452
coef exp(coef) se(coef) z Pr(>|z|)
rxLev -0.010789 0.989269 0.111379 -0.097 0.92283
rxLev+5FU -0.375746 0.686776 0.119575 -3.142 0.00168 **
age 0.007366 1.007393 0.004051 1.818 0.06901 .
surg 0.243871 1.276179 0.103202 2.363 0.01813 *
obstruct 0.283165 1.327324 0.116138 2.438 0.01476 *
differentiationpoor 0.373783 1.453222 0.121599 3.074 0.00211 **
differentiationwell 0.069022 1.071459 0.165690 0.417 0.67699
node4 0.929854 2.534140 0.098507 9.440 < 2e-16 ***
local_spreadmuscle -0.995556 0.369518 0.252106 -3.949 7.85e-05 ***
local_spreadserosa -0.501169 0.605822 0.194154 -2.581 0.00984 **
local_spreadsubmucosa -1.322018 0.266597 0.536289 -2.465 0.01370 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
exp(coef) exp(-coef) lower .95 upper .95
rxLev 0.9893 1.0108 0.79526 1.2306
rxLev+5FU 0.6868 1.4561 0.54329 0.8682
age 1.0074 0.9927 0.99943 1.0154
surg 1.2762 0.7836 1.04248 1.5623
obstruct 1.3273 0.7534 1.05711 1.6666
differentiationpoor 1.4532 0.6881 1.14506 1.8443
differentiationwell 1.0715 0.9333 0.77436 1.4826
node4 2.5341 0.3946 2.08921 3.0738
local_spreadmuscle 0.3695 2.7062 0.22545 0.6057
local_spreadserosa 0.6058 1.6506 0.41408 0.8864
local_spreadsubmucosa 0.2666 3.7510 0.09319 0.7627
Concordance= 0.672 (se = 0.012 )
Likelihood ratio test= 145.1 on 11 df, p=<2e-16
Wald test = 149.1 on 11 df, p=<2e-16
Score (logrank) test = 159.7 on 11 df, p=<2e-16anova(m3)Analysis of Deviance Table
Cox model: response is Surv(time, status)
Terms added sequentially (first to last)
loglik Chisq Df Pr(>|Chi|)
NULL -2930.2
rx -2924.1 12.1478 2 0.0023022 **
age -2923.9 0.3443 1 0.5573434
surg -2921.7 4.5404 1 0.0331029 *
obstruct -2919.3 4.7876 1 0.0286648 *
differentiation -2911.3 15.9492 2 0.0003441 ***
node4 -2867.2 88.1915 1 < 2.2e-16 ***
local_spread -2857.6 19.1615 3 0.0002532 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1cox_summary <- tidy(m3)
cox_summary# A tibble: 11 × 5
term estimate std.error statistic p.value
<chr> <dbl> <dbl> <dbl> <dbl>
1 rxLev -0.0108 0.111 -0.0969 9.23e- 1
2 rxLev+5FU -0.376 0.120 -3.14 1.68e- 3
3 age 0.00737 0.00405 1.82 6.90e- 2
4 surg 0.244 0.103 2.36 1.81e- 2
5 obstruct 0.283 0.116 2.44 1.48e- 2
6 differentiationpoor 0.374 0.122 3.07 2.11e- 3
7 differentiationwell 0.0690 0.166 0.417 6.77e- 1
8 node4 0.930 0.0985 9.44 3.74e-21
9 local_spreadmuscle -0.996 0.252 -3.95 7.85e- 5
10 local_spreadserosa -0.501 0.194 -2.58 9.84e- 3
11 local_spreadsubmucosa -1.32 0.536 -2.47 1.37e- 2Test whether proportional hazard assumptions are met for multipvariable moel predictors
m3 <- coxph(Surv(time, status) ~ rx + age + surg + obstruct +
differentiation + node4 + local_spread, data = df)
summary(m3)Call:
coxph(formula = Surv(time, status) ~ rx + age + surg + obstruct +
differentiation + node4 + local_spread, data = df)
n= 929, number of events= 452
coef exp(coef) se(coef) z Pr(>|z|)
rxLev -0.010789 0.989269 0.111379 -0.097 0.92283
rxLev+5FU -0.375746 0.686776 0.119575 -3.142 0.00168 **
age 0.007366 1.007393 0.004051 1.818 0.06901 .
surg 0.243871 1.276179 0.103202 2.363 0.01813 *
obstruct 0.283165 1.327324 0.116138 2.438 0.01476 *
differentiationpoor 0.373783 1.453222 0.121599 3.074 0.00211 **
differentiationwell 0.069022 1.071459 0.165690 0.417 0.67699
node4 0.929854 2.534140 0.098507 9.440 < 2e-16 ***
local_spreadmuscle -0.995556 0.369518 0.252106 -3.949 7.85e-05 ***
local_spreadserosa -0.501169 0.605822 0.194154 -2.581 0.00984 **
local_spreadsubmucosa -1.322018 0.266597 0.536289 -2.465 0.01370 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
exp(coef) exp(-coef) lower .95 upper .95
rxLev 0.9893 1.0108 0.79526 1.2306
rxLev+5FU 0.6868 1.4561 0.54329 0.8682
age 1.0074 0.9927 0.99943 1.0154
surg 1.2762 0.7836 1.04248 1.5623
obstruct 1.3273 0.7534 1.05711 1.6666
differentiationpoor 1.4532 0.6881 1.14506 1.8443
differentiationwell 1.0715 0.9333 0.77436 1.4826
node4 2.5341 0.3946 2.08921 3.0738
local_spreadmuscle 0.3695 2.7062 0.22545 0.6057
local_spreadserosa 0.6058 1.6506 0.41408 0.8864
local_spreadsubmucosa 0.2666 3.7510 0.09319 0.7627
Concordance= 0.672 (se = 0.012 )
Likelihood ratio test= 145.1 on 11 df, p=<2e-16
Wald test = 149.1 on 11 df, p=<2e-16
Score (logrank) test = 159.7 on 11 df, p=<2e-16anova(m3)Analysis of Deviance Table
Cox model: response is Surv(time, status)
Terms added sequentially (first to last)
loglik Chisq Df Pr(>|Chi|)
NULL -2930.2
rx -2924.1 12.1478 2 0.0023022 **
age -2923.9 0.3443 1 0.5573434
surg -2921.7 4.5404 1 0.0331029 *
obstruct -2919.3 4.7876 1 0.0286648 *
differentiation -2911.3 15.9492 2 0.0003441 ***
node4 -2867.2 88.1915 1 < 2.2e-16 ***
local_spread -2857.6 19.1615 3 0.0002532 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1cox.zph(m3) # final model with stepwise variable selection chisq df p
rx 2.3352 2 0.31111
age 0.5487 1 0.45885
surg 0.0197 1 0.88827
obstruct 6.1477 1 0.01316
differentiation 17.4588 2 0.00016
node4 5.6618 1 0.01734
local_spread 7.0826 3 0.06931
GLOBAL 37.5250 11 9.4e-05zph_test <- cox.zph(m3)
print(zph_test) chisq df p
rx 2.3352 2 0.31111
age 0.5487 1 0.45885
surg 0.0197 1 0.88827
obstruct 6.1477 1 0.01316
differentiation 17.4588 2 0.00016
node4 5.6618 1 0.01734
local_spread 7.0826 3 0.06931
GLOBAL 37.5250 11 9.4e-05# plot the Schoenfeld residuals
plot(zph_test)
Insight: Differentiation, node4 and obstruct variables did not meet proportional hazards assumption.
Stratify model by variables violating roportional hazard assumption
m4 <- coxph(Surv(time, status) ~ rx + age + surg + strata(obstruct) + strata(differentiation) + node4 +
local_spread, data = df)
summary(m4)Call:
coxph(formula = Surv(time, status) ~ rx + age + surg + strata(obstruct) +
strata(differentiation) + node4 + local_spread, data = df)
n= 929, number of events= 452
coef exp(coef) se(coef) z Pr(>|z|)
rxLev -0.019015 0.981164 0.111678 -0.170 0.86480
rxLev+5FU -0.349150 0.705288 0.119460 -2.923 0.00347 **
age 0.008636 1.008674 0.004074 2.120 0.03400 *
surg 0.258591 1.295104 0.103252 2.504 0.01226 *
node4 0.917350 2.502650 0.099125 9.255 < 2e-16 ***
local_spreadmuscle -1.067692 0.343801 0.252445 -4.229 2.34e-05 ***
local_spreadserosa -0.552546 0.575483 0.194432 -2.842 0.00449 **
local_spreadsubmucosa -1.445837 0.235549 0.536904 -2.693 0.00708 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
exp(coef) exp(-coef) lower .95 upper .95
rxLev 0.9812 1.0192 0.78828 1.2212
rxLev+5FU 0.7053 1.4179 0.55806 0.8914
age 1.0087 0.9914 1.00065 1.0168
surg 1.2951 0.7721 1.05783 1.5856
node4 2.5027 0.3996 2.06075 3.0393
local_spreadmuscle 0.3438 2.9087 0.20962 0.5639
local_spreadserosa 0.5755 1.7377 0.39313 0.8424
local_spreadsubmucosa 0.2355 4.2454 0.08224 0.6747
Concordance= 0.674 (se = 0.015 )
Likelihood ratio test= 122.9 on 8 df, p=<2e-16
Wald test = 126.4 on 8 df, p=<2e-16
Score (logrank) test = 133.9 on 8 df, p=<2e-16summary(m4)Call:
coxph(formula = Surv(time, status) ~ rx + age + surg + strata(obstruct) +
strata(differentiation) + node4 + local_spread, data = df)
n= 929, number of events= 452
coef exp(coef) se(coef) z Pr(>|z|)
rxLev -0.019015 0.981164 0.111678 -0.170 0.86480
rxLev+5FU -0.349150 0.705288 0.119460 -2.923 0.00347 **
age 0.008636 1.008674 0.004074 2.120 0.03400 *
surg 0.258591 1.295104 0.103252 2.504 0.01226 *
node4 0.917350 2.502650 0.099125 9.255 < 2e-16 ***
local_spreadmuscle -1.067692 0.343801 0.252445 -4.229 2.34e-05 ***
local_spreadserosa -0.552546 0.575483 0.194432 -2.842 0.00449 **
local_spreadsubmucosa -1.445837 0.235549 0.536904 -2.693 0.00708 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
exp(coef) exp(-coef) lower .95 upper .95
rxLev 0.9812 1.0192 0.78828 1.2212
rxLev+5FU 0.7053 1.4179 0.55806 0.8914
age 1.0087 0.9914 1.00065 1.0168
surg 1.2951 0.7721 1.05783 1.5856
node4 2.5027 0.3996 2.06075 3.0393
local_spreadmuscle 0.3438 2.9087 0.20962 0.5639
local_spreadserosa 0.5755 1.7377 0.39313 0.8424
local_spreadsubmucosa 0.2355 4.2454 0.08224 0.6747
Concordance= 0.674 (se = 0.015 )
Likelihood ratio test= 122.9 on 8 df, p=<2e-16
Wald test = 126.4 on 8 df, p=<2e-16
Score (logrank) test = 133.9 on 8 df, p=<2e-16cox.zph(m4) # final model with stratification by proportional hazard violating variables chisq df p
rx 2.0007 2 0.368
age 0.6704 1 0.413
surg 0.0142 1 0.905
node4 4.2882 1 0.038
local_spread 5.2976 3 0.151
GLOBAL 12.4113 8 0.134zph_test <- cox.zph(m4)
print(zph_test) chisq df p
rx 2.0007 2 0.368
age 0.6704 1 0.413
surg 0.0142 1 0.905
node4 4.2882 1 0.038
local_spread 5.2976 3 0.151
GLOBAL 12.4113 8 0.134# plot the Schoenfeld residuals
plot(zph_test)
Insight: After model stratification by obstruct and differentiation, the proportional hazard assumption is met.
Perform 5-fold cross-validation for model evaluation with proper stratification of boot samples
To include proper stratification in all 5 samples, including fold assignment steps
library(survival)
library(boot) # for bootstrapping
library(survcomp) # to calculate c-index
library(caret)
# set a seed
set.seed(1234)
# Cox model
cox_model <- coxph(Surv(time, status) ~ rx + age + surg + strata(obstruct) + strata(differentiation) + node4 + local_spread, data = df)
# calculate the original c-index
c_index_original <- concordance.index(predict(cox_model), surv.time = df$time, surv.event = df$status)$c.index
cat("original c-index:", c_index_original, "\n")original c-index: 0.6545844 # create a function for calculating c-index in each fold
cox_cindex <- function(train_data, test_data) {
fit <- coxph(Surv(time, status) ~ rx + age + surg + strata(obstruct) + strata(differentiation) + node4 + local_spread, data = train_data)
predictions <- predict(fit, newdata = test_data)
c_index <- concordance.index(predictions, surv.time = test_data$time, surv.event = test_data$status)$c.index
return(c_index)
}
# perform 5-fold cross-validation with stratification
K <- 5
folds <- createFolds(c(df$status, df$differentiation, df$rx), k = K, list = TRUE, returnTrain = TRUE)
cv_c_indices <- sapply(folds, function(train_indices) {
train_data <- df[train_indices, ]
test_data <- df[-train_indices, ]
cox_cindex(train_data, test_data) # cox_cindex is defined to fit the cox model on training data and calculate the c-index for test data in each fold.
}) # createFolds from the caret package to perform stratified 5-fold cross-validation, ensuring balance across strata
# cross-validation c-indices
cat("cross-validated c-Indices for each fold:", cv_c_indices, "\n")cross-validated c-Indices for each fold: 0.6852433 0.6165189 0.6402717 0.6272685 0.7060098 cat("mean cross-validated c-Index:", mean(cv_c_indices), "\n")mean cross-validated c-Index: 0.6550624 # plot cross-validation c-indices
plot(cv_c_indices, type = "b", xlab = "Fold", ylab = "c-index", main = "c-index across folds")
Insight: The original model c-index (0.654) and mean cross-validation c-index (0.655) is very similar, suggesting the the final model is stable and is not overfitting.
V. Results
Table 2. Univariate model: Survival after Chemotherapy for stage B/C Colon Cancer
| Treatment | Coefficient | Hazard ratio | 95% CI_upper | 95% CI_lower | P-value |
|---|---|---|---|---|---|
| Amisole (Lev) | -0.027 | 0.974 | 0.784 | 1.209 | 0.809 |
| Amisole + 5-FU | -0.372 | 0.690 | 0.546 | 0.870 | 0.002 |
Table 3. Multivariate model: Survival after Chemotherapy for stage B/C Colon Cancer
| Treatment | Coefficient | Hazard ratio | 95% CI_upper | 95% CI_lower | P-value |
|---|---|---|---|---|---|
| Amisole (Lev) | -0.011 | 0.989 | 0.795 | 1.231 | 0.923 |
| Amisole + 5-FU | -0.376 | 0.687 | 0.543 | 0.868 | 0.002 |
| Age | 0.007 | 1.007 | 0.999 | 1.015 | 0.069 |
| Surge | 0.244 | 1.276 | 1.042 | 1.562 | 0.018 |
| Obstruction of colon | 0.283 | 1.327 | 1.057 | 1.667 | 0.015 |
| Differentiat ion_poor | 0.374 | 1.453 | 1.145 | 1.844 | 0.002 |
| Differentiat ion_well | 0.069 | 1.072 | 0.774 | 1.483 | 0.677 |
| More than 4 nodes (+) | 0.930 | 2.534 | 2.089 | 3.074 | 3.75 x 10-21 |
| Local spread_muscle | -0.996 | 0.370 | 0.225 | 0.606 | 7.85 x 10-5 |
| Local spread_serosa | -0.501 | 0.606 | 0.414 | 0.886 | 0.010 |
| Local spread_submucosa | -1.322 | 0.267 | 0.093 | 0.763 | 0.014 |
References
[1]
Terry M. Therneau and Patricia M. Grambsch, Modeling survival data: Extending the Cox model. New York: Springer, 2000.
[2]
T. M. Therneau, A package for survival analysis in r. 2024. Available: https://CRAN.R-project.org/package=survival
[3]
R Core Team, R: A language and environment for statistical computing. Vienna, Austria: R Foundation for Statistical Computing, 2023. Available: https://www.R-project.org/
[4]
A. Kassambara, M. Kosinski, and P. Biecek, Survminer: Drawing survival curves using ’ggplot2’. 2021. Available: https://CRAN.R-project.org/package=survminer
[5]
W. N. Venables and B. D. Ripley, Modern applied statistics with s, Fourth. New York: Springer, 2002. Available: https://www.stats.ox.ac.uk/pub/MASS4/
[6]
A. C. Davison and D. V. Hinkley, Bootstrap methods and their applications. Cambridge: Cambridge University Press, 1997. Available: doi:10.1017/CBO9780511802843
[7]
Angelo Canty and B. D. Ripley, Boot: Bootstrap r (s-plus) functions. 2024.